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I. INTRODUCTION 

Propagation of the cosmic rays in the turbulent mag¬ 
netic fields can proceed in different regimes depending on 
the scales under consideration. On small scales, when the 
particles move coherently, their propagation is ballistic. 
This usually occurs close to the source, just after the par¬ 
ticles escape the sites of their acceleration. With time, 
the multiple stochastic scattering in turbulent magnetic 
fields leads to the isotropization of directions of cosmic 
rays. The complete isotropization implies that the prop¬ 
agation proceeds in the diffusion regime. 

The limiting cases of small-angle and isotropic particle 
distributions allow solutions of the problem of particle 
propagation on small and large scales, respectively Una. 
The small-angle approximation fails when the deflection 
of particles becomes large, typically larger than one ra¬ 
dian. The solution of the diffusion equation, in addition 
to its inability to be applied to small spacial scales, faces 
the so-called problem of superluminal propagation [4|. 
In this regard, the apparent requirement r 2 /D > r/c im¬ 
plies that the diffusion works when r > D/c, where D is 
the diffusion coefficient. 

In the small-angle approximation, the evolution of an¬ 
gular distribution has a diffusive behaviour mm- How¬ 
ever, since the pitch angle changes within the limited 
interval, — 1 < p < 1 (/i is the cosine of the pitch angle), 
the mean square displacement of the pitch angle deviates 
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from diffusive behaviour as the average deflection angle 
grows. For the isotropic turbulence, the moments of the 
pitch-angle distribution on large timescales have an expo¬ 
nential behaviour [5] . This means that the isotropization 
becomes fast after the characteristic time which is deter¬ 
mined by the pitch-angle Fokkcr-Planck coefficient D 
On the other hand, a slower isotropization can happen if 
the turbulent magnetic field has a slab geometry. 

In this work we use the method of moments. It allows 
us to eliminate the angular dependence at the expense 
of introduction of an isotropization function that deter¬ 
mines the dynamics of isotropization. The final results 
depend on the form of the isotropization function. How¬ 
ever, as we demonstrate below, for a reasonable choice of 
the form of the function, the results remain quite stable. 

To avoid of the problem of the superluminal motion, 
Aloisio and Berezinsky Hj have introduced the so-called 
Jiittner function, which describes the evolution of the 
cosmic ray density. Although this function is obtained 
phenomenologically from the formal similarity between 
the diffusion propagator and the Maxwellian distribution, 
it gives correct results in the limiting cases of diffusion 
and the ballistic regime. Below we will show that the in¬ 
tegrated over time Jiittner function is close our stationary 
solution, which proceeds from Boltzmann equation. 

In many cases, especially for the problems related to 
the radiation of cosmic rays, it is necessary to know not 
only the cosmic-ray density but also their angular distri¬ 
bution. Indeed, the part of radiation emitted by particles 
with strongly anisotropic distribution, can have a strong 
impact on the morphology and spectrum of radiation, or 
even simply missed by the observer. 
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To demonstrate the importance of the angular distribu¬ 
tion for calculations of the apparent gamma-ray morphol¬ 
ogy, we calculated the gamma-ray intensity maps of the 
regions surrounding the cosmic ray accelerator. The most 
distinct features can be seen for the clumps of matter 
close to the source. The anisotropy changes significantly 
the radiation spectrum and leads to fast decrease of in¬ 
tensity at high energies. This results in a suppression or 
a disappearance of radiation from nearby clouds located 
away from the line of sight towards to the cosmic-ray 
source. Moreover, even in the case of homogeneous mat¬ 
ter distribution, the effects of anisotropy and the transi¬ 
tion from ballistic to diffusion regime, play an important 
role in the formation of gamma-ray morphology. 

The article is organised as follows. In Section [II] we de¬ 
scribe the formalism based on the Boltzmann equation, 
and obtain the stationary distribution function, which is 
valid for all scales from the ballistic to diffusion regime. 
In Section |III[ this distribution function is used for cal¬ 
culations of the gamma-ray morphology of nearby clouds 
irradiated by cosmic rays. The conclusions are summa¬ 
rized in Section m 


II. ANALYTICAL DESCRIPTION 
A. Method of moments 


Let us consider the evolution of the distribution func¬ 
tion in the case of multiple stochastic scatterings. Here 
we do not specify the mechanism of the scattering, and 
describe the processes only by the generic probability of 
the particle to be scattered. The evolution of the dis¬ 
tribution function f(t,r,n) is determined by the Boltz¬ 
mann transport equation 


df Of 

at + "S? = st/ 


8{t)S(r) 
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(1) 


where the speed of the light (propagation speed of ultra- 
relativistic particles) is set c = 1. Here n is the unit vec¬ 
tor in the direction of propagation. It is assumed that the 
particles are produced by an instant spherically symmet¬ 
ric point-like source described by the Dirac delta func¬ 
tions S(t)6(r). We consider only elastic collisions which 
are described by the collision integral 


Stf = j f(n')w(n' -A n)dVL' — J f(ri)w{n -A n')d£t', 


( 2 ) 

where w(n' -A n) is the probability of scattering of the 
particle from the initial direction along n' to the final 
direction along n per unit time. In the case of isotropic 
medium, the probability w depends only on the angle be¬ 
tween the initial and final directions. For compactness of 
presentation, the dependence of the distribution function 
on time and coordinates in Eq. 0 is omitted. 

The solution of Eq. 0 can be found in the small-angle 
approximation which is valid for initial moments of time. 


Below we will show that proceeding from this equation 
one can obtain also the equation for diffusion of parti¬ 
cles, and derive its solution which is valid for large time 
intervals. 

We are interested, first of all, in the transition between 
these two solutions. For these purpose it is useful to sim¬ 
plify the problem and consider instead of the distribution 
function its moments. Applying successively the integral 
operators f dfl, f dfln a , ... , f dQn a ...n UJ to the Boltz¬ 
mann equation, one can obtain the equations for the mo¬ 
ments. We restrict ourselves to the first two moments: 
the density g = f fdtt and the current j = f nfdfl. 
They are governed by the following equations: 


dg dj a 

dt dr a 
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Here r is the scattering time, which is the inverse of the 
transport cross section cr tr = 1/r, where 

a tr = J (1 — nn')w(n' —■> n)dVL. (4) 

In the derivation of equations in Eq. ([3]) it has been taken 
into account that J n a dQ = 0, / Stf dfl = 0. 

The density and the current depend on the higher mo¬ 
ments of the distribution function. To close the system 
of equations given by Eq. |3|, we should define the form 
of the isotropization tensor, 

, \ / n a npfdn 

{nan0} = Jfdn ’ (5) 

based on the following physical arguments. 

In the spherically symmetric case, the radial direc¬ 
tion is the only preferential direction, therefore, the 
isotropization tensor should have the following structure 


(' n a np) = AS a p + Bp a pp , (6) 

where p = r/r is the radial direction. The standard 
procedure for the determination of the coefficients A and 
B, which consists in the consequent multiplication by 
tensors 5 a p and p a pp, leads to the equations 

1 = 3H + B (7) 

(( np ) 2 ) = A + B, 


from where we find 

A= 1 1 M) i g _ 3((rep) 2 ) - 1 (g) 

We assume that the tensor depends only on the coordi¬ 
nates. Then, it is convenient to introduce the isotropiza¬ 
tion function <f>(r) = B which would allow us to write the 
isotropization tensor in the form 

(n a np) = (1 - 0( r ))^ + <f)(r)p a pp. (9) 
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The tensor consists of the unidirectional, p a pp, and 
the isotropic, 8 a p/3, parts. At r = 0, when n = p, 
we have (n a np}(0) = p a pp, whereas at the infinity 
all directions of n are distributed isotropically, and 
(n a np)(oo) = 5 a p/3. Thus, the isotropization function 
should satisfy the boundary conditions 

<f)(0) = 1 and 4>(oo) = 0. (10) 


In the spherically symmetric case, the density and the 
current are expressed as g = g(t 7 r) and j = j(t 7 r)p. 
Then the substitution of the isotropization tensor in the 
form given by Eq. © into Eq. © results in 


dG dJ_ 

dt dr 

dJ dG 

dt dr 


8(t)S(r) 
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( 11 ) 


where we have introduced the functions G = r 2 g and 
J = r 2 j. 

At small distances and times, the current changes fast, 
i.e. ^ 3> y. Therefore one can neglect the term in the 
right-hand side of the second equation in Eq. (11). The 


condition </>(0) = 1 leads to the cancellation of the third 
term in the left-hand side of the same equation. Thus, 


the equations in Eq. (Ill are reduced to 


dG dJ__ 8{t)5[r) 
dt ^ dr 47t 

dJ dG _ 

~dt + Ih ~ ’ 


( 12 ) 


which can be rewritten in the form of the wave equation 
for G 


d 2 G d 2 G 
~dt? ~ ~d^’ 


(13) 


with the boundary condition G(t = +0, r) = <5(r)/47r. 
The solution of this equation in terms of the density g = 
G/r 2 is 


9^ r ) = ^2 6 (r-t), (14) 


which describes the behaviour of density in the ballistic 
regime. 

In the opposite case , i.e. for large distances and times, 
it is more convenient to proceed from the initial system of 
equations given by Eq. The current changes slowly 
on scale of the scattering time r, i.e. ^ A This allows 
us to neglect the derivative ^ in the second equation. 
Taking into account that <j> = 0 we find 



j = ~^g, 


(15) 


which can be rewritten in the form of the diffusion equa¬ 
tion for g 


dg _ t 

dt 3 Aff ' 


(16) 


The comparison with the conventional form of diffusion 
equation, gives the well known relation 


D = 


Cl c 

T’ 


(17) 


where l c = ct is the scattering length. 


B. Stationary case 


The system of equations derived in the previous sec¬ 
tion describes the particle motion over the entire pro¬ 
cess of propagation, including the ballistic and diffusion 
modes, as well as the transition stage between these two 
regimes. They have a simple stationary solution. Indeed, 
the integration over the entire time cancels out the time 
derivatives and leads to the following system of ordinary 
differential equations 


dJ 

dr 

dG 

dr 


6{r) 
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where J = J(r) and G = G(r). The first of these equa¬ 
tions gives J = 0(r)/47r, where 0(r) is the Heaviside 
step function. The substitution of the current J to the 
second equation results in the ordinary differential equa¬ 
tion of the first order. It’s solution can be presented in 
the form 



The limits for the integrals are chosen from the condition 
that the density should vanish in the infinity. 

To find the final expression, one should choose a suit¬ 
able form of the isotropization function <j>(r). The bound¬ 
ary condition 4irg(r)r 2 —> 1 at r —> 0, which has not been 
used yet, gives a relation between <p(r) and r. To take 
into account this relation, </> should have one free param¬ 
eter v. In the case of isotropic medium, one expects an 
exponential rate of isotropization [5]. Therefore, one of 
the possible expressions for isotropization function is 

</>(r) = e~ r/v . 


( 21 ) 












However, this function does not allow representation of 
Eq. (19) in quadratures. Another expression 
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1+r/v 


( 22 ) 


is less physically motivated, but allows us to obtain a 
simple analytical solution. Indeed, the substitution of 
Eq. (22) into Eq. (191 results in 


d(r) = 


1 3 (u + r) 

4 7 T7’2 t 


(23) 


The boundary condition 47rg(r)r 2 —> 1 at r — > 0 gives 
v = t/3. The comparison with Eq. ( [l7| gives v = D /c. 
Then, the solution can be written in terms of diffusion 
coefficient D as 


5(r) 47 t (r 2 t\d) ’ ^ 

which is just the sum of the solutions in the limiting cases 
of ballistic and diffusion regimes. Introducing a dimen¬ 
sionless parameter x = rc/D, which is the radial distance 
in the units of D /c, one can rewrite this expression in the 
form 


9(r) = 


0-+x) 

Airr 2 


(25) 


The result is surprisingly simple. However, we 
should note that it is obtained for a specific form of 
isotropization function. The results for other forms of 
isotropization function can be numerically calculated us¬ 
ing Eq. (191. Fig. [l] shows that the densities obtained 
for <fi = and = e _r /" differ less than 30%, 

although it should be noticed that in the first case the 
asymptotic behaviour occurs later. For comparison, we 
present also the density provided by the generalized 
Jiittner function integrated over time. This function was 
proposed phenomenologically in ref. |4] for description of 
evolution of the cosmic ray flux. It is seen from Fig.[j]that 
the Jiittner function gives a result which is quite close to 
the solution with an exponential form of the isotropiza¬ 
tion function. 

The propagation of cosmic rays is characterized not 
only by the density, but also by the angular distribution 
of particles. Such an information is contained in the mo¬ 
ments of the angular distribution. Indeed, if p = np is 
the cosine of the angle between the particle and the radial 
direction, then the first moment of the angular distribu¬ 
tion is 


(p) = 


3P 
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Using the solution given by Eq. (251, we have 


(p) 


1 

1 + x' 


(26) 


(27) 



Figure 1: Upper panel: Density G = 4nr 2 g(r) for different 
models of isotropization as function of x = cr/D. Gi(x ) and 
Gi{x) correspond to (j >i = 1 and (p 2 = e~ rv isotropization 
functions. Gj(x) corresponds to the integrated over time gen¬ 
eralized Jiittner function proposed by jTE]. The dashed lines 
represents the asymptotes G(x) = 1 and G{x) = x. Lower 
panel: The ratios of the functions: G 2 /G 1 and Gj/Gi. 


The second moment is just the projection of the 
isotropization tensor on the radial direction 


(p 2 ) = {n a np)p a p p = 1 + . (28) 

The moments (p) and (p 2 ) can be used to construct a 
function that has properties of the exact angular distri¬ 
bution. In this regard, the simplest function is 

M (p) = ^ exp ’ ( 29 ) 

where the normalization function has the form 

Z(aO = !(l-e- fl /*). (30) 


M(p) is the distribution of the cosine of pitch angles p 
in the range [—1,1], This distribution gives the moments 
that are in a good agreement, at small and large dis¬ 
tances, with moments given by Eqs. (271 and (281. For 
both moments the difference from Eqs. (27) and (281 does 
not exceed 35%. 

One can show that in the small-angle limit the distri¬ 
bution given by Eq. (29) describes the diffusion in angle. 
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Indeed, in this limit the transport cross section becomes 

ft i n 0\. 1A6 2 Dg 

a tr « y (1 — (1 — y))w(n -A n)dO = -<-^-) = ~Y’ 

(31) 

where Dg is the diffusion coefficient in angular space. 
Taking into account Eq. ( fl7| and the relation x = rc/D, 
we obtain 

M(ri ~ exp (- = exp (- -AAj . (32) 


The combination of the density and the angular distri¬ 
butions results in the stationary distribution function 




Q ( 1 c\ 1 

47TC V f 2 tD ) 27 iZ ^ 


3Z?(1 — n) 


rc 


(33) 

which describes the evolution of propagation, from bal¬ 
listic to diffusion, where Q is the source function (the 
production rate) of cosmic rays in the source. 


C. Diffusion coefficient 


leads to the dependence of the diffusion coefficient on 
the turbulence spectrum. The quasi-linear theory pre¬ 
dicts l c ~ [7], where a is the power-law index of the 

turbulence spectrum. This dependence can be obtained 
as follows. The scattering frequency of the particle by 
waves is expressed as |8j 



S’kk \ 

B 2 /^) \k=k reB ’ 


(37) 


where $k is the spectrum of turbulence normalised as 
f $kdk = B 2 /8n, k res = f l/vfj, is the resonance wave 
vector, is the girofrequency, and fj, is the cosine of the 
pitch angle. The diffusion coefficient is related to the 
scattering frequency as mm 



(38) 


Taking turbulence spectrum in the form $k ~ l/k a 
with minimum wave vector ko = 27r/A, and substitut¬ 
ing Eq. (37) into Eq. (381, we obtain 


In the derived distribution functions the diffusion co¬ 
efficient is assumed to be constant in space, but it can 
have an arbitrary energy dependence. The latter is 
determined by the relation between the Larmor radius 
Rl = E/eB and the correlation length of the turbulence 
A. If Rl A, the particle is only slightly deflected on the 
correlation length. The random walk of uncorrelated de¬ 
flections results in the scattering length l c ~ R 2 L By 
combining Eq. Eq. pTj ) and the relation r = l/<r tr , 
we obtain 


2 c 2 

3 2V 


(34) 


The coefficient of diffusion in angle, Dg, can be written 

as |2 


Dg 


(a - 1)(/? + 1) cA 
Ia/3 W L ' 


(35) 


where a and /3 are the power-law indices of turbulence 
for wave vectors k > ko = 2 tt/\ and k < ko, respectively. 
The value of (3 is poorly known; here we will assume 
0 = 1 . 

For the Kolmogorov spectrum of turbulence (a = 5/3), 
calculations of the diffusion coefficient in the regime 
Rl > A give 



(36) 


Note that the energy dependence in this regime is the 
same for any other spectrum of turbulence. 

At R^ <C A, particles are only scattered by MHD 
waves with the length equal to their gyroradius. This 


n= _2_ 1 Rl 

(2 - a)(4 - a) [a - 1) 7 r(2 7 r)«- 1 \ A 


2—a 


Thus for a = 5/3 we have 


D = 


27 


7 7t(27t) 2 / 3 


Ac 


f) 


1/3 


Rl < A. 


(39) 


(40) 


We note that the value of the numerical prefactors in 
Eqs. (36) and (401 depend slightly on the on turbulence 


spectrum. The parameter A in these expressions corre¬ 
sponds to the largest scale of inertial range of turbulence 
and can be associated with the correlation length. 

The simulations of the particle propagation in the 
isotropic and purely turbulent magnetic field performed 
in ref. [5] show that the diffusion coefficient can be pre¬ 
sented in the following form 



1 


Rl\ 1/3 
(2tt) 2 /3 V A ) 



+wHe ■ ( 4 i) 


In the next section, this convenient expression for the 
diffusion coefficient will be used for calculations of the 
gamma-ray emission. 


III. GAMMA-RAY EMISSION OF THE REGION 
SURROUNDING COSMIC-RAY ACCELERATOR 

High energy gamma-rays carry unique information not 
only about the accelerators of cosmic rays (electrons, pro¬ 
tons and nuclei), but also allow us to trace these particle 
after they leave the sites of their acceleration. In the 
interstellar medium, this is realized through interactions 
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of cosmic rays with the so-called giant molecular clouds. 
These dense gas regions illuminated by cosmic rays pro¬ 
vide a target for proton-proton interaction and radiate 
gamma-rays. Thus they can serve as unique "barome¬ 
ters" for measurements of the pressure (energy density) 
of cosmic rays at different locations relative to the accel¬ 
erator. 

The massive clouds located in the vicinity of the accel¬ 
erator dramatically increase the chances of tracing the 
run-away particles through the secondary gamma-rays. 
For example, for a young supernova remnant (as an ac¬ 
celerator of cosmic rays) at a distance of 1 kpc, a gas 
cloud of mass of order 10 4 Mq can emit very high energy 
gamma-rays at a level detectable by current instruments 
if the cloud is located within 100 pc from the supernova 
remnant m■ Before being fully diffused way and inte¬ 
grated into the "sea" of the galactic cosmic rays, they 
produce gamma-rays the spectrum of which could essen¬ 
tially differ from both the gamma-ray spectrum of the 
accelerator itself and the spectrum of the diffuse galac¬ 
tic gamma-ray emission. In the case of propagation in a 
"nominal" diffusion regime, the formation of gamma-ray 
energy spectra has been discussed in ref. m- However, 
closer to the accelerator, the propagation of cosmic rays 
may have a more complex character which would be re¬ 
flected in the spectra of secondary gamma-rays. 

The morphology of the gamma-ray emission is deter¬ 
mined by the interplay between cosmic ray and mat¬ 
ter distributions. While the distribution of the matter 
could be arbitrary and random, the density of the cos¬ 
mic rays decreases with distance from the accelerator, 
which means the most intensive radiation should arrive 
from the dense regions located close to the accelerator. 

However, the angular distribution of cosmic rays may 
lead to a significant deviation from such a simple picture 
m- Indeed, due to the relativistic character of proton- 
proton interactions, gamma rays are emitted along the 
direction of the momentum of the incident proton. It 
means that only the protons directed towards the ob¬ 
server give a contribution to the detectable gamma-ray 
flux. Close to the accelerator, the angular distribution of 
cosmic rays, which propagate ballistically along the ra¬ 
dial direction, is strictly anisotropic. Thus the gamma 
radiation can arrive from the direction towards accelera¬ 
tor only in the case of presence of significant amount of 
gas along the line of sight. The nearby clouds located 
not on the line of sight could be invisible because their 
radiation is directed not towards the observer. As the 
distribution of cosmic rays becomes more isotropic, the 
apparent intensity of the radiation increases with the dis¬ 
tance from the accelerator. 

On the other hand, if a gas cloud on the line of sight is 
located sufficiently close to the accelerator, we could see a 
bright source which would coincide with the accelerator. 
This can be misinterpreted as a prolific production of 
gamma-rays inside the accelerator, although in reality the 
accelerator could be a very inefficient gamma-ray emitter. 

Thus, the consideration of transition from the ballistic 


to diffusion regime may result in two consequences re¬ 
lated to the density and angular distributions of cosmic 
rays. The first one is the lack of radiation from regions 
close to the accelerator even in the case of presence of 
massive nearby clouds (but located away from the line 
of sight). In contrary, we may detect very bright and fo¬ 
cused gamma ray image if a dense cloud in proximity of 
the accelerator would appear on the line of sight. 

In the following calculations we consider a pro¬ 
ton accelerator of power L cr = 10 3 'erg/s located 

at a distance d = 1 kpc from the observer. The 
energy spectrum of protons is taken in the form 
J cr (E) = E~ 2 exp (—.E/10 15 eV). The stationary dis¬ 
tribution of cosmic rays around the accelerator is de¬ 
scribed by Eq. ( [33] ) with the diffusion coefficient given by 
Eq. (411. We assume that particles propagate through a 
turbulent magnetic field with the Kolmogorov spectrum 
of turbulence. We keep the coherence length of turbu¬ 
lence A as a free parameter. 

In calculations, the magnetic field is taken at the 
level of B = 10 ” 4 G. In the galactic disk, the typi¬ 
cal value of the diffusion coefficient is D ss 10 2 S cm 2 /s 
at 1 GeV H3]. To match this value we as¬ 
sume A = 10 4 pc and A = 10 5 pc for the diffu¬ 

sion coefficients, which below we refer to as small 
and large, respectively. The smallness of the ratio 
Rl/ A ~ 1.1 x 10 _ 9 Ei 2 AJ 1 , where E 12 = E/ 10 12 eV is 
the cosmic ray energy and A 4 = A/10 4 pc, allows us to 
neglect the high-energy non-resonant part of the diffu¬ 
sion coefficient. The estimate based on Eq. (41) gives 
D zz 9 x lO^E^A ^ 3 cm 2 /s. Correspondingly, the tran¬ 
sition from ballistic to diffusion regime occurs at a charac¬ 
teristic distance D/c « 1.0 -E^A ^ 3 pc from the cosmic- 
ray source. 

Below we consider two cases: (i) a homogeneous cloud 
surrounding the accelerator, (ii) a group of clouds lo¬ 
cated near the accelerator. We consider different den¬ 
sities of the background gas (between the clouds), and 
different values of the diffusion coefficient of cosmic rays. 
For calculation of gamma-ray production in pp collisions 
the parametrisation of ref. m has been used. The re¬ 
sults present the intensity of gamma-radiation integrated 
along the line of sight through the gamma-ray production 
region. 

In the first example, we consider a homogeneous cloud 
with a radius R = 10 pc and density n p = 100cm” 3 . At 
these parameters the mass of the cloud is M c ; ss 10 4 Mq. 
The accelerator is placed in the centre of the cloud. The 
intensity maps of gamma rays at three different energies, 
E = lO 10 eV, E = 10 11 eV, and E = 10 12 eV, are shown 
in Fig. [2] Calculations are performed for the case of 
slow diffusion. It is assumed that the density inside the 
accelerator is very small, thus the gamma-radiation of 
the accelerator itself can be ignored. 

Fig. [ 2 ] shows that with the increase of energy the radi¬ 
ation becomes more concentrated on the position of the 
accelerator. The change of the gamma-ray morphology 
is explained by the energy dependence of the cosmic-ray 
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propagation. With the increase of energy, the transition 
from ballistic to diffusion regime occurs at larger dis¬ 
tances from the accelerator. The faster decrease of the 
cosmic-ray density in the ballistic regime results in the 
faster decrease of gamma-ray intensity at high energies. 
It is seen in Fig. [3]in which the radial profiles of intensity 
at different energies are shown. The left and right pan¬ 
els of Fig. [3] present the results calculated for the small 
and large diffusion coefficients, respectively. For the fast 
diffusion, the transition from the steep to the flat part of 
the profile occurs at large distances from the accelerator. 

In Fig. [3] we show the radial profiles of gamma-ray 
intensities for two different regimes of diffusion. To de¬ 
scribe the behaviour of these curves, we also show the 
local power-law indices a = d In P(r)/dln(r) at different 
energies. One can see that at small distances all profiles 
approach to the same inclination with a ~ —1.3. At 
small energies the flat part starts earlier and intensity 
decreases slower. For example, at energy E = 10 9 eV the 
flattest part of the profile corresponds toaw —0.3. With 
an increase of energy, the flat part becomes steeper. At 
very high energies the transition to the flat part might 
not happen at all. 


The angular distribution of cosmic rays has a strong 
impact on the characteristics of the secondary gamma- 
radiation, in particular on the radial profile of the 
gamma-ray intensity (see Appendix 0. In the case of 
strictly radial distribution of particles (i.e. for the the 
ballistic regime of propagation), the gamma-ray source 
will be detected as a point like object, independent of 
the linear size of the gamma-ray production region, L. 
For the “nominal” diffusion regime, the angular size of 
the gamma-ray source is determined by the ratio L/R. 
The transition regime introduces non-neglible corrections 
to the formation of the overall image of the gamma-ray 
source, therefore should be treated thoroughly. The com¬ 
parison of the results in Fig. [3] and Fig. A.l (which does 
not take into account the angular distribution of protons) 
shows that, in general, the shapes of the radial profiles 
in these figures are similar, but there is also a significant 
difference. In particular, if one assumes that angular dis¬ 
tribution of particles is isotropic just after their escape 
from the source,the power-law index of the slope of the 
(projected) gamma-ray profile would be a = — 1. This 
is in contrast to a ss —1.3, which is expected if we cor¬ 
rectly treat the angular distribution of particles closer 
to the source. The sharp angular distribution of cosmic 
rays close to the accelerator leads to considerable loss of 
emission, and therefore to a steeper gamma-ray intensity 
profile. 

The spectral energy distributions (SED) of gamma 
rays at different distances from the source of cosmic rays 
are shown in Fig. [4] The left and right panels show the 
results corresponding to the small and large diffusion co¬ 
efficients, respectively. It is seen that the increase of the 
diffusion coefficient leads to harder gamma-ray spectra. 
In the case of homogeneous distribution of the cosmic-ray 
density and for the power-law energy spectrum of protons 


with an index a = 2, gamma-rays have an almost flat 
SED. The flat part slightly deviates from the cosmic ray 
spectrum, namely it contains an intrinsic hardening due 
to the increase of inelastic cross section m■ For the ho¬ 
mogeneous distribution of cosmic rays, the flat part of the 
gamma-ray spectrum can be approximated by power-law 
with an photon index T = 1.94. The SED in the direction 
to accelerator becomes even harder with the power-law 
photon index F = 1.86. This can be explained by the 
fact that high-energy protons preserve longer their ra¬ 
dial direction and, consequently, produce higher energy 
radiation towards observer. 

Because of diffusion of cosmic rays, in the directions 
far from direction to the accelerator, the density of the 
low-energy protons decreases slower. At large distances 
from the accelerator, the SED of gamma-rays is close to 
E~ 5 (slightly harder because of the pp interaction cross- 
section), where 6 characterizes the energy-dependence of 
the diffusion coefficient, D(E) = DqE s . Because of the 
assumed dependence of the diffusion coefficient with <5 = 
1/3, the SED at large distances in Fig. [4] follows oc E 1 / 3 
behaviour. 

In more realistic scenarios, the surroundings of the 
cosmic-ray source could be inhomogeneous, i.e. may con¬ 
sist of clumps of matter. To study the general features 
of radiation of such an environment, we use a simpli¬ 
fied gas distribution-template consisting of four identical 
equally separated clouds surrounded by a homogeneous 
low-density gas (background): 



where n p 0 = 10 3 cm -3 , rq and Wi are the coordinates of 
the centres of these clouds and their widths, respectively; 
Xb g is the level of the background relative to the maxi¬ 
mum density in the centres of the clouds. The width of 
each cloud is Wj = 1 pc which corresponds to the mass 
M c i « 140 Mq. Separation between clouds is 5 pc. We 
consider the case without background and with the back¬ 
ground of the level of Xb g = 10 -2 . The zero level of the 
background allows us to eliminate the radiation from the 
direction towards the accelerator. For illustrative pur¬ 
poses, the source of cosmic rays is located at the left 
border of each map in Fig. [8] and Fig. [9j It is assumed 
that the density inside the accelerator is very low, thus 
its own gamma-radiation can be neglected. 

The intensity maps for the case without background 
gas calculated for the slow and fast diffusion of cosmic 
rays, are shown in Figs. [5] and [6] respectively. The corre¬ 
sponding intensity profiles are presented in the left and 
right panels of Fig. [7} For slow diffusion, the gamma-ray 
intensity of clouds decreases with the distance from the 
proton accelerator. The clouds are located along the line 
perpendicular to the line of sight. At high energies in the 
case of fast diffusion a part of the radiation from the clos¬ 
est clouds is "lost". Therefore, one can see an interesting 
effect when despite the decrease of the cosmic ray density 



Table I: Power-law index a for the fit to the gamma-ray in¬ 
tensity profiles for positions of the maximum radia tion (the 
centres of the clouds) shown in the right panel of Fig. |10| (large 
diffusion coefficient). 


-E 7 , eV 

10 9 

10 10 

to 11 

10 12 

10 13 

10 14 

a 

-0.85 

-0.73 

-0.49 

-0.21 

0.53 

1.14 


with distance, highest energy gamma-rays are seen from 
the furthest rather than closest clouds. 

The results corresponding to the homogeneous gas 
background are shown in Figs. [8] and [9] The intensity 
profiles shown in Fig. [10] are smoother compared to the 
relevant curves in the case of background absence. The 
radiation towards the accelerator appears in Fig. [10] be¬ 
cause of the presence of the background gas on the line 
of sight. The results in Fig. [T0| show that in the case of 
fast diffusion the gamma-ray intensity is reduced; at very 
high energies it disappears at all from the closest cloud. 

To describe quantitatively how fast the gamma-ray in¬ 
tensity decreases with distance, the intensity at the posi¬ 
tion of the maximum (in the centres of the clouds) shown 
in Fig. [TO] has been fitted with power law. The calcu¬ 
lated power-law indices of the fits for different energies 
are presented in Table [IJ It is seen that with increase 
of energy the profile becomes flatter. The intensity de¬ 
creases slower in the case of large diffusion coefficient. 
Moreover, for energies E — 10 13 and E = 10 14 the inten¬ 
sity increases with distance for the farthest two clouds. In 
all cases the intensity in maximum points changes with 
distance slower than 1/p, where p is the projected dis¬ 
tance from the source. 


The spectral energy distributions of gamma rays in 
the direction to the centres of the clouds are shown in 
Figs. 12 and 11 The clouds are numbered in the order of 
their distances to the accelerator. The energy spectra of 
gamma-rays from the clouds are steeper than in the case 
of the homogeneous cloud surrounding the accelerator 
(see Fig. [l|. Obviously, this is explained by anisotropic 
distribution of cosmic rays closer to the accelerator. The 
picture becomes smoother if the clouds are embedded in 
a homogeneous background gas. 


IV. CONCLUSION 


One of the major issues in the general problem of iden¬ 
tification of sources of galactic and extragalactic cosmic 
rays, is the character of their propagation through the 
turbulent magnetic fields outside the accelerators. De¬ 
pending on the distance to the source, the level of the 
magnetic turbulence of the ambient medium, as well as 
the energy of particles, their propagation can proceed in 
the ballistic or in the diffusion regimes. While the diffu¬ 
sion of cosmic rays has been comprehensively studies in 
the literature, the description of propagation in the in¬ 


termediate stage, i.e. at the transition from the ballistic 
to the diffusive regime, is a problem of greater complex¬ 
ity regarding the exact analytical solutions. To study 
the dynamics of this transition, we used a simplified ap¬ 
proach. Based on the treatment of moments of the Boltz¬ 
mann equation, we derived the system of equations for 
the time evolution of the distribution function of cosmic 
rays. Written in the form of Eq. (Ill, it describes not 
only the ballistic and diffusion regimes of propagation 
as limiting cases, but also the transition stage between 
these two regimes. This system of equations allows a 
simple stationary solution in the form of Eq. (19). 

The key feature of this approach is the proper 
choice of a specific function which describes the pitch- 
angle isotropization. Within the chosen approach, the 
isotropization function cannot be strictly determined. 
Nevertheless it can be chosen and introduced in a self- 
consistent manner based on reasonable physics argu¬ 
ments. In this paper, two forms of the isotropization 
function has been considered: (1) (f> = e~ r ^ u and (2) 
(f> = 1/(1 + r/is). The exponential form seems to be bet¬ 
ter justified given that at the presence of isotropic turbu¬ 
lence the pitch angle distribution moments behave expo¬ 
nentially [5]. On the other hand the form <f> = 1/(1 +r/v) 
provides a simple representation of the stationary solu¬ 
tion given by Eq. (23 ). The results obtained for both 
forms of the isotropization function are in agreement 
within better than 30% accuracy. They agree well also 
with the integrated Jiittner function proposed in ref. [3]. 

The angular, energy and radial distributions of cosmic 
rays outside the accelerator, in the general case of their 
propagation, including the transition stage between the 
ballistic and diffusion regimes, is described by a surpris¬ 
ingly simple function given by Eq. (33 ). 

The anisotropy of angular distribution of cosmic rays 
before they enter the diffusion regime of propagation, 
leads to a partial or complete "loss" of secondary gamma- 
rays by the observer. The numerical calculations of sec¬ 
ondary gamma-rays performed for both the homogeneous 
and clumpy (consisting of dense clouds) distributions of 
gas in the vicinity of the accelerator, demonstrate the 
strong impact of the effects of particle propagation in 
the pre-diffusion regime on the apparent gamma-ray mor¬ 
phology and the energy spectrum. Therefore, the de¬ 
tailed studies of spectral and morphological features of 
high energy gamma-ray emission outside the detected (or 
potential) cosmic ray accelerators can tell us about the 
propagation character of cosmic rays after they leave the 
sites of their acceleration. 


Appendix A: Intensity profiles 

The gamma-ray intensity profile corresponds to the ra¬ 
diation integrated along the line of sight and considered 
as a function of projected distance from the source. Let 
us assume that the angular distribution of cosmic-ray 
protons is isotropic. In this case the production rate of 
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Figure A.l: The gamma-ray intensity profiles calculated 
for a spherical homogeneous cloud surrounding the source 
of protons. The curves correspond to three values of the 
parameter X. The upper panel: the intensity, lower panel 
- the slop of the profile. 



Figure A.2: The same as in Fig. |A.1| but for a strip 
region. 


gamma rays in the direction of the observer is propor¬ 
tional to the density of protons. Since the cosmic-ray 
density changes according to Eq. (24), the form of inten¬ 
sity profile is 




(Al) 


where z axis is directed along the line of sight. 

If the region emitting gamma-rays is a spherical 
cloud of radius R with the centre at the position 
of the particle accelerator, then the integration limits 
Zi ,2 = ±\/R 2 3 4 - p 2 give 


J-sph — — I 


/ arctan 


- 1 


R 


4- Xarccosh 





(A2) 


where p is the projected distance from the source, p = 
p/R, and X = Rc/D. If the emitting region is a strip 
with width 2 R, and the source of cosmic rays is in the 


middle of the strip, then the integration limits Z\^ = ±R 
give 


2 (arctan 1 

- _ _ I _ v_ 

str Ry p 


A'arcsinh 



(A3) 


The functions I sp hR/2 and I s t r R /2 are presented 
in Fig. B and Fig. B respectively. The lower 
panels show the local slops of the curves defined as 
a = d\nF{rj)/d\n(rj), where F(p) = IR/2. 

It is seen from Figs. |A.l] and |A.2] that at small distances 
the intensity changes as 1/p. For X = 1 this behaviour 
approximately retains at all distances. This curve cor¬ 
responds to the case when the diffusion regime has not 
been reached at the distance R. Other curves reveal the 
transition to the diffusion regime becoming flatter with 
increase of distance. The integration over smaller region 
along the line of sight in the case of a spherical cloud 
makes the slope of the profile steeper compared to the 
case of strip region. 
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Intensity map, erg/cm 2 s sr 



Figure 2: The intensity maps of gamma-ray emission at different energies. The spherical cloud with homogeneous density 
distribution is irradiated by the cosmic-ray source located in its centre. The gas density inside the accelerator is assumed very 
low, so the contribution of the accelerator to the gamma-ray emission is negligible. The maps are produced for the case of small 
diffusion coefficient (for details, see the text). For the distance to the source d = 1 kpc, the region of ~ 1° x 1° corresponds to 
the area ~ 20 x 20 pc 2 . 



Figure 3: Upper panel: radial intensity profiles for different energies of gamma-rays in the case of homogeneous cloud 
surrounding the cloud. Lower panel: the power-law index of the corresponding intensity profiles in the upper panel. The left 
and right panels present the results calculated for the case of small and large diffusion coefficients, respectively. The profiles 
for energies E = Id 10 eV, E = 10 11 eV, and E = 10 12 eV in the left panel correspond to the intensity maps shown in Fig. [ 2 ] 
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Figure 4: Energy spectra of gamma rays at different distances from the cosmic-ray source in the case of homogeneous cloud. 
The results are given for the scenarios with slow (left panel) and fast (right panel) proton diffusion. 
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Figure 5: The intensity maps of gamma-ray emission from the group of clouds (without background) at various energies for 
the case of low diffusion coefficient. The cosmic-ray source is located in the centre of the left side. 
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Figure 6: The same as in Fig. [5] for the case of high diffusion coefficient. 
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Figure 7: The radial intensity profiles for various energies in the case of group of clouds for the low (left panel) and high (right 
panel) diffusion coefficient. The profiles for energies E = 10 9 eV, E = 10 11 eV, and E = 10 13 eV correspond to the intensity 
maps given in Fig. [5] for the left panel and Fig. [6] for the right panel. 
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Figure 8: The intensity maps of gamma-ray emission from the group of clouds and homogeneous background at various 
energies for the case of low diffusion coefficient. The cosmic-ray source is located in the centre of the left side. 
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Figure 9: The same as in Fig. [8] for the case of high diffusion coefficient. 
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Figure 10: The radial intensity profiles for various energies in the case of group of clouds surrounded by homogeneous 

background for the low (left panel) and high (right panel) diffusion coefficient. The profiles for energies E = 10 9 eV, E = 10 11 eV, 
and E = 10 13 eV correspond to the intensity maps given in Fig. [8] for the left panel and Fig. [ 9 ] for the right panel. 



Figure 11: Energy spectra of gamma rays in the direction to the centres of the clouds in the case of background absence for 
low (left panel) and high (right panel) diffusion coefficient. The clouds are numbered in the order of the distance from the 
source. 
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Figure 12: Energy spectra of gamma rays in the direction to the centres of the clouds in the case of homogeneous background 
for low (left panel) and high (right panel) diffusion coefficient. The clouds are numbered in the order of the distance from the 
source. 




